import numpy as np
from scipy.sparse import linalg, csr_matrix
import config
import math

class TwoFluidSolver:
    def __init__(self, subchannel):
        """初始化两流体模型求解器
        
        Args:
            subchannel: SubChannel对象，包含子通道的几何和物性信息
        """
        self.subchannel = subchannel
        
        # 网格参数
        self.n_axial = subchannel.n_axial
        self.n_channels = subchannel.n_channels
        self.dx = config.GEOMETRY_PARAMS['channel_length'] / (subchannel.n_axial - 1)
        self.dy = config.GEOMETRY_PARAMS['channel_width']
        
        # 时间步长控制
        self.dt = config.NUMERICAL_PARAMS['time_step']
        self.dt_min = config.NUMERICAL_PARAMS['min_time_step']
        self.dt_max = config.NUMERICAL_PARAMS['max_time_step']
        self.dt_increase = config.NUMERICAL_PARAMS['time_step_increase_factor']
        self.dt_decrease = config.NUMERICAL_PARAMS['time_step_decrease_factor']
        self.total_time = config.NUMERICAL_PARAMS['total_time']
        self.current_time = 0.0
        
        # 收敛控制
        self.max_iterations = config.NUMERICAL_PARAMS['max_iterations']
        self.convergence_tolerance = config.NUMERICAL_PARAMS['convergence_tolerance']
        
        # 松弛因子
        self.relaxation_factor = config.NUMERICAL_PARAMS['relaxation_factor']
        self.lateral_flow_under_relaxation = config.NUMERICAL_PARAMS['lateral_flow_under_relaxation']
        self.pressure_under_relaxation = config.NUMERICAL_PARAMS['pressure_under_relaxation']
        self.void_fraction_under_relaxation = config.NUMERICAL_PARAMS['void_fraction_under_relaxation']
        
        # 初始化存储上一步解的数组
        self._initialize_previous_solution()
        
    def _initialize_previous_solution(self):
        """初始化存储上一步解的数组"""
        shape = (self.n_axial, self.n_channels)
        self._previous_pressure = np.zeros(shape)
        self._previous_void_fraction = np.zeros(shape)
        self._previous_vel_liquid = np.zeros(shape)
        self._previous_vel_vapor = np.zeros(shape)
        self._previous_temp_liquid = np.zeros(shape)
        self._previous_temp_vapor = np.zeros(shape)
        self._previous_lateral_velocity = np.zeros(shape)
        
    def solve_transient(self):
        """求解瞬态问题"""
        n_timesteps = int(self.total_time / self.dt)
        
        try:
            for step in range(n_timesteps):
                self.current_time = (step + 1) * self.dt
                print(f"\n开始计算时间步 {step+1}/{n_timesteps}, 当前时间: {self.current_time:.3f}s")
                print(f"当前时间步长: {self.dt:.2e}s")
                
                # 保存上一时间步的解
                self._store_previous_solution()
                
                # 求解当前时间步
                converged = self._solve_current_timestep()
                
                if not converged:
                    # 减小时间步长并重试
                    self.dt = max(self.dt * self.dt_decrease, self.dt_min)
                    print(f"未收敛，减小时间步长至: {self.dt:.2e}s")
                    if self.dt == self.dt_min:
                        raise RuntimeError("时间步长已达到最小值仍未收敛")
                    # 恢复上一步的解
                    self._restore_previous_solution()
                    step -= 1
                    continue
                    
                # 如果收敛，可以适当增加时间步长
                self.dt = min(self.dt * self.dt_increase, self.dt_max)
                
                # 更新边界条件
                self._update_boundary_conditions()
                
                # 定期输出结果
                if (step + 1) % config.OUTPUT_PARAMS['print_frequency'] == 0:
                    self._print_results()
                    
        except Exception as e:
            print(f"计算过程中出错: {str(e)}")
            raise
            
    def _solve_current_timestep(self):
        """求解当前时间步
        
        Returns:
            bool: 是否收敛
        """
        for iter in range(self.max_iterations):
            # 求解动量方程
            self._solve_momentum_equations()
            
            # 求解质量方程
            self._solve_mass_equations()
            
            # 求解能量方程
            self._solve_energy_equations()
            
            # 更新本构关系
            self._update_constitutive_relations()
            
            # 检查收敛性
            if self._check_convergence():
                print(f"迭代收敛，迭代次数: {iter+1}")
                return True
                
        print(f"警告：达到最大迭代次数 {self.max_iterations} 仍未收敛")
        return False
        
    def _check_convergence(self):
        """检查收敛性
        
        Returns:
            bool: 是否收敛
        """
        # 检查压力场的收敛性
        dp_max = np.max(np.abs(self.subchannel.pressure - self._previous_pressure))
        if dp_max > self.convergence_tolerance:
            return False
            
        # 检查速度场的收敛性
        du_l_max = np.max(np.abs(self.subchannel.vel_liquid - self._previous_vel_liquid))
        du_v_max = np.max(np.abs(self.subchannel.vel_vapor - self._previous_vel_vapor))
        if max(du_l_max, du_v_max) > self.convergence_tolerance:
            return False
            
        # 检查空泡分数的收敛性
        da_max = np.max(np.abs(self.subchannel.void_fraction - self._previous_void_fraction))
        if da_max > self.convergence_tolerance:
            return False
            
        # 检查温度场的收敛性
        dt_l_max = np.max(np.abs(self.subchannel.temp_liquid - self._previous_temp_liquid))
        dt_v_max = np.max(np.abs(self.subchannel.temp_vapor - self._previous_temp_vapor))
        if max(dt_l_max, dt_v_max) > self.convergence_tolerance:
            return False
            
        return True
        
    def _store_previous_solution(self):
        """存储上一时间步的解，用于时间导数项的计算"""
        # 存储速度场
        self._previous_vel_liquid = self.subchannel.vel_liquid.copy()
        self._previous_vel_vapor = self.subchannel.vel_vapor.copy()
        self._previous_lateral_velocity = self.subchannel.lateral_velocity.copy()  # 修正：lateral_flow -> lateral_velocity
        
        # 存储压力场
        self._previous_pressure = self.subchannel.pressure.copy()
        
        # 存储空泡分数
        self._previous_void_fraction = self.subchannel.void_fraction.copy()
        
        # 存储温度场
        self._previous_temp_liquid = self.subchannel.temp_liquid.copy()
        self._previous_temp_vapor = self.subchannel.temp_vapor.copy()
        self._previous_temp_wall = self.subchannel.temp_wall.copy()
        
    def _restore_previous_solution(self):
        """恢复上一次迭代的解"""
        self.subchannel.pressure = self._previous_pressure.copy()
        self.subchannel.void_fraction = self._previous_void_fraction.copy()
        self.subchannel.vel_liquid = self._previous_vel_liquid.copy()
        self.subchannel.vel_vapor = self._previous_vel_vapor.copy()
        self.subchannel.temp_liquid = self._previous_temp_liquid.copy()
        self.subchannel.temp_vapor = self._previous_temp_vapor.copy()
        self.subchannel.lateral_velocity = self._previous_lateral_velocity.copy()
        
    def _print_results(self):
        """打印计算结果"""
        print("\n当前计算结果:")
        print(f"最大压力: {np.max(self.subchannel.pressure):.2e} Pa")
        print(f"最大空泡分数: {np.max(self.subchannel.void_fraction):.3f}")
        print(f"最大液相速度: {np.max(self.subchannel.vel_liquid):.2f} m/s")
        print(f"最大气相速度: {np.max(self.subchannel.vel_vapor):.2f} m/s")
        print(f"最大液相温度: {np.max(self.subchannel.temp_liquid):.2f} K")
        print(f"最大气相温度: {np.max(self.subchannel.temp_vapor):.2f} K")
        
    def _update_boundary_conditions(self):
        """更新边界条件"""
        # 入口边界条件
        self.subchannel.pressure[0,:] = config.INITIAL_CONDITIONS['inlet_pressure']
        self.subchannel.temp_liquid[0,:] = config.INITIAL_CONDITIONS['inlet_temperature']
        self.subchannel.temp_vapor[0,:] = config.INITIAL_CONDITIONS['inlet_temperature']
        self.subchannel.void_fraction[0,:] = config.INITIAL_CONDITIONS['inlet_void_fraction']
        
        # 出口边界条件（零梯度）
        self.subchannel.pressure[-1,:] = self.subchannel.pressure[-2,:]
        self.subchannel.temp_liquid[-1,:] = self.subchannel.temp_liquid[-2,:]
        self.subchannel.temp_vapor[-1,:] = self.subchannel.temp_vapor[-2,:]
        self.subchannel.void_fraction[-1,:] = self.subchannel.void_fraction[-2,:]
        
    def _solve_momentum_equations(self):
        """求解动量方程"""
        # 轴向动量方程
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                self._solve_axial_momentum('l', i, j)  # 液相
                self._solve_axial_momentum('v', i, j)  # 气相
                
        # 横向动量方程
        for i in range(self.n_axial):
            for j in range(self.n_channels - 1):
                self._solve_lateral_momentum('l', i, j)  # 液相
                self._solve_lateral_momentum('v', i, j)  # 气相
                
    def _solve_axial_momentum(self, k, i, j):
        """求解轴向动量方程 (eq.2)"""
        if k == 'l':
            alpha = 1 - self.subchannel.void_fraction[i,j]
            rho = self.subchannel.rho_liquid[i,j]
            u = self.subchannel.vel_liquid[i,j]
            u_old = self.subchannel.vel_liquid[i,j]  # 使用当前值作为初始值
        else:
            alpha = self.subchannel.void_fraction[i,j]
            rho = self.subchannel.rho_vapor[i,j]
            u = self.subchannel.vel_vapor[i,j]
            u_old = self.subchannel.vel_vapor[i,j]  # 使用当前值作为初始值
            
        # 防止alpha为0或1时的除零错误
        alpha = max(1e-6, min(1-1e-6, alpha))
        rho = max(rho, 1e-6)  # 防止密度为0
            
        # 如果存在上一时间步的值，则使用它
        if hasattr(self, '_previous_vel_liquid') and k == 'l':
            u_old = self._previous_vel_liquid[i,j]
        elif hasattr(self, '_previous_vel_vapor') and k == 'v':
            u_old = self._previous_vel_vapor[i,j]
            
        # 间导数项
        d_time = alpha * rho * (u - u_old) / self.dt
        
        # 对流项（一阶迎风格式）
        if i > 0:
            if u > 0:
                u_up = self.subchannel.vel_liquid[i-1,j] if k == 'l' else self.subchannel.vel_vapor[i-1,j]
            else:
                u_up = u
        else:
            u_up = u
            
        if i < self.subchannel.n_axial - 1:
            if u > 0:
                u_down = u
            else:
                u_down = self.subchannel.vel_liquid[i+1,j] if k == 'l' else self.subchannel.vel_vapor[i+1,j]
        else:
            u_down = u
            
        d_conv = alpha * rho * u * (u_down - u_up) / self.dx
        
        # 压力梯度
        if i < self.subchannel.n_axial - 1:
            dp_dx = (self.subchannel.pressure[i+1,j] - self.subchannel.pressure[i,j]) / self.dx
        else:
            dp_dx = (self.subchannel.pressure[i,j] - self.subchannel.pressure[i-1,j]) / self.dx
            
        # 重力项
        gravity = alpha * rho * config.PHYSICAL_PARAMS['gravity']
        
        # 壁面摩擦力
        f_wall = self._calculate_wall_friction(k, i, j)
        
        # 界面剪切力
        f_interface = self._calculate_interfacial_friction(i, j)
        
        # 相变动量
        gamma = self._calculate_phase_change(k, i, j)
        m_phase = u * gamma
        
        # 计算新的速度（添加数值稳定性保护）
        denominator = alpha * rho
        if abs(denominator) < 1e-6:
            acceleration = 0.0
        else:
            acceleration = (-d_conv - alpha * dp_dx - gravity - f_wall + f_interface + m_phase) / denominator
            
        # 限制加速度的大小以提高数值稳定性
        max_acceleration = 1000.0  # 最大允许加速度
        acceleration = max(min(acceleration, max_acceleration), -max_acceleration)
        
        # 更新速度
        if k == 'l':
            self.subchannel.vel_liquid[i,j] = u + acceleration * self.dt
        else:
            self.subchannel.vel_vapor[i,j] = u + acceleration * self.dt
            
    def _solve_lateral_momentum(self, k, i, j):
        """求解横向动量方程 (eq.3)
        
        Args:
            k: 相态 ('l' 或 'v')
            i, j: 网格位置
        """
        if j >= self.subchannel.n_channels - 1:
            return  # 最后一个通道没有横向流动
            
        if k == 'l':
            alpha = 1 - self.subchannel.void_fraction[i,j]
            rho = self.subchannel.rho_liquid[i,j]
            v = self.subchannel.lateral_velocity[i,j] / (rho * alpha) if alpha > 0 else 0
            v_old = self._get_previous_value('lateral_velocity', i, j) / (rho * alpha) if alpha > 0 else 0
        else:
            alpha = self.subchannel.void_fraction[i,j]
            rho = self.subchannel.rho_vapor[i,j]
            v = self.subchannel.lateral_velocity[i,j] / (rho * alpha) if alpha > 0 else 0
            v_old = self._get_previous_value('lateral_velocity', i, j) / (rho * alpha) if alpha > 0 else 0
            
        # 时间导数项
        d_time = alpha * rho * (v - v_old) / self.dt
        
        # 横向压力梯度
        dp_dy = (self.subchannel.pressure[i,j+1] - self.subchannel.pressure[i,j]) / self.dy
        
        # 横向壁面阻力
        f_wall_y = self._calculate_wall_friction_lateral(k, i, j)
        
        # 横向界面剪切力
        f_interface_y = self._calculate_interfacial_friction_lateral(i, j)
        
        # 相变动量
        gamma = self._calculate_phase_change(k, i, j)
        m_phase_y = v * gamma
        
        # 湍流混合动量
        t_momentum = self._calculate_turbulent_momentum(k, i, j)
        
        # 计算新的横向速度
        acceleration = (-alpha * dp_dy - f_wall_y + f_interface_y - m_phase_y + t_momentum) / (alpha * rho)
        v_new = v + acceleration * self.dt
        
        # 更新横向流量
        if k == 'l':
            self.subchannel.lateral_velocity[i,j] = v_new * rho * alpha
        else:
            self.subchannel.lateral_velocity[i,j] += v_new * rho * alpha  # 累加气相贡献
        
    def _solve_mass_equations(self):
        """求解质量守恒方程"""
        for i in range(self.n_axial):
            for j in range(self.n_channels):
                # 计算相变率
                gamma = self._calculate_phase_change('l', i, j)
                
                # 计算横向质量流量
                w_lateral = 0.0
                if j < self.n_channels - 1:
                    w_lateral += self.subchannel.lateral_velocity[i,j]
                if j > 0:
                    w_lateral -= self.subchannel.lateral_velocity[i,j-1]
                    
                # 计算轴向质量流量
                if i > 0:
                    w_axial_in_l = (1 - self.subchannel.void_fraction[i-1,j]) * \
                                 self.subchannel.rho_liquid[i-1,j] * \
                                 self.subchannel.vel_liquid[i-1,j]
                    w_axial_in_v = self.subchannel.void_fraction[i-1,j] * \
                                 self.subchannel.rho_vapor[i-1,j] * \
                                 self.subchannel.vel_vapor[i-1,j]
                else:
                    # 口界条件
                    w_axial_in_l = config.INITIAL_CONDITIONS['inlet_mass_flux'] * \
                                 (1 - config.INITIAL_CONDITIONS['inlet_void_fraction'])
                    w_axial_in_v = config.INITIAL_CONDITIONS['inlet_mass_flux'] * \
                                 config.INITIAL_CONDITIONS['inlet_void_fraction']
                    
                w_axial_out_l = (1 - self.subchannel.void_fraction[i,j]) * \
                              self.subchannel.rho_liquid[i,j] * \
                              self.subchannel.vel_liquid[i,j]
                w_axial_out_v = self.subchannel.void_fraction[i,j] * \
                              self.subchannel.rho_vapor[i,j] * \
                              self.subchannel.vel_vapor[i,j]
                
                # 更新空泡分数
                A = config.GEOMETRY_PARAMS['flow_area']
                dt = self.dt
                
                d_alpha = (w_axial_in_v - w_axial_out_v + gamma) * dt / (A * self.subchannel.rho_vapor[i,j])
                self.subchannel.void_fraction[i,j] += self.void_fraction_under_relaxation * d_alpha
                
                # 限制空泡分数在[0,1]范围内
                self.subchannel.void_fraction[i,j] = max(0.0, min(1.0, self.subchannel.void_fraction[i,j]))
                
    def _solve_energy_equations(self):
        """求解能量方程 (eq.3)"""
        # 从配置文件获取几何参数
        A = config.GEOMETRY_PARAMS['flow_area']
        P_heated = config.GEOMETRY_PARAMS['heated_perimeter']
        
        for i in range(self.subchannel.n_axial):
            for j in range(self.subchannel.n_channels):
                # 获取当前网格的物性参数
                alpha_l = max(1e-6, 1 - self.subchannel.void_fraction[i,j])  # 防止除零
                alpha_v = max(1e-6, self.subchannel.void_fraction[i,j])
                rho_l = max(1e-6, self.subchannel.rho_liquid[i,j])
                rho_v = max(1e-6, self.subchannel.rho_vapor[i,j])
                cp_l = self.subchannel.cp_liquid[i,j]
                cp_v = self.subchannel.cp_vapor[i,j]
                
                # 计算传热系数
                h_l = self._calculate_heat_transfer_coefficient('l', i, j)
                h_v = self._calculate_heat_transfer_coefficient('v', i, j)
                
                # 计算壁面传热
                T_w = self.subchannel.temp_wall[i,j]
                T_l = self.subchannel.temp_liquid[i,j]
                T_v = self.subchannel.temp_vapor[i,j]
                
                # 使用配置文件中的加热周长
                q_l = h_l * P_heated * (T_w - T_l)
                q_v = h_v * P_heated * (T_w - T_v)
                
                # 计算相变传热（分别计算液相和气相的相变）
                gamma_l = self._calculate_phase_change('l', i, j)  # 液相相变率
                gamma_v = self._calculate_phase_change('v', i, j)  # 气相相变率
                h_fg = self._calculate_latent_heat(self.subchannel.pressure[i,j])  # 潜热
                q_phase = (gamma_l + gamma_v) * h_fg  # 总相变传热
                
                # 计算温度变化
                # 添加数值保护，防止除零
                denominator_l = alpha_l * rho_l * cp_l * A
                denominator_v = alpha_v * rho_v * cp_v * A
                
                if denominator_l > 1e-6:
                    dT_l = (q_l - q_phase) * self.dt / denominator_l
                else:
                    dT_l = 0.0
                    
                if denominator_v > 1e-6:
                    dT_v = (q_v + q_phase) * self.dt / denominator_v
                else:
                    dT_v = 0.0
                
                # 更新温度
                self.subchannel.temp_liquid[i,j] += dT_l
                self.subchannel.temp_vapor[i,j] += dT_v
                
    def _update_constitutive_relations(self):
        """更新本构关系"""
        for i in range(self.subchannel.n_axial):
            for j in range(self.subchannel.n_channels):
                # 更新饱和温度和压力关系
                P = self.subchannel.pressure[i,j]
                T_sat = self._calculate_saturation_temperature(P)
                
                # 更新液相物性
                T_l = self.subchannel.temp_liquid[i,j]  # 修正：liquid_temp -> temp_liquid
                self.subchannel.rho_liquid[i,j] = self._calculate_liquid_density(T_l, P)
                self.subchannel.mu_liquid[i,j] = self._calculate_liquid_viscosity(T_l)
                self.subchannel.k_liquid[i,j] = self._calculate_liquid_conductivity(T_l)
                self.subchannel.cp_liquid[i,j] = self._calculate_liquid_specific_heat(T_l)
                
                # 更新气相物性
                T_v = self.subchannel.temp_vapor[i,j]  # 修正：vapor_temp -> temp_vapor
                self.subchannel.rho_vapor[i,j] = self._calculate_vapor_density(T_v, P)
                self.subchannel.mu_vapor[i,j] = self._calculate_vapor_viscosity(T_v)
                self.subchannel.k_vapor[i,j] = self._calculate_vapor_conductivity(T_v)
                self.subchannel.cp_vapor[i,j] = self._calculate_vapor_specific_heat(T_v)
                
                # 更新潜热
                self.subchannel.latent_heat[i,j] = self._calculate_latent_heat(T_sat)
                
                # 更新表面张力
                self.subchannel.surface_tension[i,j] = self._calculate_surface_tension(T_sat)
                
    def _calculate_liquid_density(self, T, P):
        """计算液相密度"""
        # 简化模型：线性关系
        rho_ref = 1000.0  # 参考密度 (kg/m³)
        beta = 2e-4  # 热膨胀系数 (1/K)
        return max(1e-6, rho_ref * (1 - beta * (T - 273.15)))
        
    def _calculate_vapor_density(self, T, P):
        """计算气相密度"""
        # 理想气体方程
        R = 461.5  # 水蒸气气体常数 (J/kg·K)
        return max(1e-6, P / (R * T))
        
    def _calculate_liquid_viscosity(self, T):
        """计算液相粘度"""
        # 简化模型：指数关系
        mu_ref = 1e-3  # 参考粘度 (Pa·s)
        E_mu = 2000  # 活化能 (J/mol)
        R = 8.314  # 气体常数 (J/mol·K)
        return max(1e-6, mu_ref * math.exp(E_mu/R * (1/T - 1/373.15)))
        
    def _calculate_vapor_viscosity(self, T):
        """计算气相粘度"""
        # 简化模型：Sutherland公式
        mu_ref = 1.8e-5  # 参考粘度 (Pa·s)
        T_ref = 373.15  # 参考温度 (K)
        S = 650  # Sutherland常数 (K)
        return max(1e-6, mu_ref * (T/T_ref)**1.5 * (T_ref + S)/(T + S))
        
    def _calculate_liquid_conductivity(self, T):
        """计算液相导热系数"""
        # 简化模型：线性关系
        k_ref = 0.6  # 参考导热系数 (W/m·K)
        dk_dT = -0.0015  # 温度系数 (W/m·K²)
        return max(1e-6, k_ref + dk_dT * (T - 373.15))
        
    def _calculate_vapor_conductivity(self, T):
        """计算气相导热系数"""
        # 简化模型：线性关系
        k_ref = 0.025  # 参考导热系数 (W/m·K)
        dk_dT = 7e-5  # 温度系数 (W/m·K²)
        return max(1e-6, k_ref + dk_dT * (T - 373.15))
        
    def _calculate_liquid_specific_heat(self, T):
        """计算液相比热容"""
        # 简化模型：常数
        return 4200.0  # J/kg·K
        
    def _calculate_vapor_specific_heat(self, T):
        """计算气相比热容"""
        # 简化模型：常数
        return 2000.0  # J/kg·K
        
    def _calculate_latent_heat(self, T_sat):
        """计算汽化潜热 (J/kg)"""
        # 使用Watson公式
        T_crit = 647.15  # 临界温度 (K)
        h_fg_ref = 2257e3  # 在1个大气压下的汽化潜热 (J/kg)
        
        # 防止数值不稳定
        T_sat = min(T_sat, T_crit - 0.1)  # 确保不会超过临界温度
        T_ratio = max(0.0, min(1.0, (T_crit - T_sat)/(T_crit - 373.15)))  # 限制在[0,1]范围内
        
        return max(1e-6, h_fg_ref * T_ratio**0.38)
                
    def _calculate_surface_tension(self, T):
        """计算表面张力"""
        # 简化模型：线性关系
        sigma_ref = 0.059  # 参考表面张力 (N/m)
        dsigma_dT = -1.8e-4  # 温度系数 (N/m·K)
        return max(1e-6, sigma_ref + dsigma_dT * (T - 373.15))
        
    def _get_previous_value(self, field, i, j):
        """获取上一时间步的值"""
        if field == 'vel_liquid':
            return self._previous_vel_liquid[i,j]
        elif field == 'vel_vapor':
            return self._previous_vel_vapor[i,j]
        elif field == 'lateral_velocity':  # 修正：lateral_flow -> lateral_velocity
            return self._previous_lateral_velocity[i,j]
        elif field == 'pressure':
            return self._previous_pressure[i,j]
        elif field == 'void_fraction':
            return self._previous_void_fraction[i,j]
        elif field == 'temp_liquid':
            return self._previous_temp_liquid[i,j]
        elif field == 'temp_vapor':
            return self._previous_temp_vapor[i,j]
        elif field == 'temp_wall':
            return self._previous_temp_wall[i,j]
        else:
            raise ValueError(f'Unknown field: {field}')
        
    def _calculate_wall_friction(self, k, i, j):
        """计算壁面摩擦力 (eq.13-18)
        
        Args:
            k: 相态 ('l' 或 'v')
            i, j: 网格位置
            
        Returns:
            壁面摩擦力
        """
        if k == 'l':
            alpha = 1 - self.subchannel.void_fraction[i,j]
            rho = self.subchannel.rho_liquid[i,j]
            u = self.subchannel.vel_liquid[i,j]
            mu = self.subchannel.mu_liquid[i,j]
        else:
            alpha = self.subchannel.void_fraction[i,j]
            rho = self.subchannel.rho_vapor[i,j]
            u = self.subchannel.vel_vapor[i,j]
            mu = self.subchannel.mu_vapor[i,j]
            
        # 防止alpha为0或1时的除零错误
        alpha = max(1e-6, min(1-1e-6, alpha))
            
        # 计算雷诺数
        d_h = config.GEOMETRY_PARAMS['hydraulic_diameter']
        Re = rho * abs(u) * d_h / mu
        
        # 计算摩擦系数
        if Re < 2300:  # 层流
            f = 64 / max(Re, 1e-6)  # 防止Re=0时除零
        else:  # 湍流
            f = 0.316 / (max(Re, 1e-6)**0.25)
            
        # 计算相流修正系数
        if k == 'l':
            phi_l2 = (1 - alpha)**(-1.75)  # 液相两相流修正系数
            f *= phi_l2
        else:
            phi_g2 = alpha**(-1.75)  # 气相两相流修正系数
            f *= phi_g2
            
        # 计算壁面摩擦力
        tau_w = 0.5 * f * rho * u * abs(u)
        
        # 考虑湿周
        P_w = config.GEOMETRY_PARAMS['wetted_perimeter']
        A_c = config.GEOMETRY_PARAMS['flow_area']
        
        return tau_w * P_w / A_c
        
    def _calculate_wall_friction_lateral(self, k, i, j):
        """计算横向壁面摩擦力
        
        Args:
            k: 相态 ('l' 或 'v')
            i, j: 网格位置
            
        Returns:
            横向壁面摩擦力
        """
        if k == 'l':
            alpha = 1 - self.subchannel.void_fraction[i,j]
            rho = self.subchannel.rho_liquid[i,j]
            v = self.subchannel.lateral_velocity[i,j] / (rho * alpha) if alpha > 0 else 0
        else:
            alpha = self.subchannel.void_fraction[i,j]
            rho = self.subchannel.rho_vapor[i,j]
            v = self.subchannel.lateral_velocity[i,j] / (rho * alpha) if alpha > 0 else 0
            
        # 使用横向流动形状损失系数
        K = config.LATERAL_FLOW_PARAMS['form_loss_coefficient']
        
        return 0.5 * K * rho * v * abs(v)
        
    def _calculate_interfacial_friction(self, i, j):
        """计算界面剪切力 (eq.19)
        
        Args:
            i, j: 网格位置
            
        Returns:
            界面剪切力
        """
        # 获取相速度
        u_l = self.subchannel.vel_liquid[i,j]
        u_v = self.subchannel.vel_vapor[i,j]
        
        # 获取液相密度
        rho_l = self.subchannel.rho_liquid[i,j]
        
        # 计算界面面积
        alpha_l = 1 - self.subchannel.void_fraction[i,j]
        alpha_v = self.subchannel.void_fraction[i,j]
        d_h = config.GEOMETRY_PARAMS['hydraulic_diameter']
        A_i = 4 * alpha_l * alpha_v / d_h
        
        # 计算界面阻力系数
        C_i = config.LATERAL_FLOW_PARAMS['crossflow_resistance']
        
        return C_i * rho_l * A_i * (u_l - u_v) * abs(u_l - u_v)
        
    def _calculate_interfacial_friction_lateral(self, i, j):
        """计算横向界面剪切力
        
        Args:
            i, j: 网格位置
            
        Returns:
            横向界面剪切力
        """
        # 获取液相和气相的横向速度
        alpha_l = 1 - self.subchannel.void_fraction[i,j]
        alpha_v = self.subchannel.void_fraction[i,j]
        rho_l = self.subchannel.rho_liquid[i,j]
        
        v_l = self.subchannel.lateral_velocity[i,j] / (rho_l * alpha_l) if alpha_l > 0 else 0
        v_v = self.subchannel.lateral_velocity[i,j] / (self.subchannel.rho_vapor[i,j] * alpha_v) if alpha_v > 0 else 0
        
        # 计算界面面积
        d_h = config.GEOMETRY_PARAMS['hydraulic_diameter']
        A_i = 4 * alpha_l * alpha_v / d_h
        
        # 计算界面阻力系数
        C_i = config.LATERAL_FLOW_PARAMS['crossflow_resistance']
        
        return C_i * rho_l * A_i * (v_l - v_v) * abs(v_l - v_v)
        
    def _calculate_turbulent_momentum(self, k, i, j):
        """计算湍流混合动量
        
        Args:
            k: 相态 ('l' 或 'v')
            i, j: 网格位置
            
        Returns:
            湍流混合动量
        """
        if j >= self.subchannel.n_channels - 1:
            return 0
            
        # 获取湍流混合系数
        beta = config.LATERAL_FLOW_PARAMS['mixing_coefficient']
        
        # 计算混合速度
        if k == 'l':
            v1 = self.subchannel.vel_liquid[i,j]
            v2 = self.subchannel.vel_liquid[i,j+1]
            rho = self.subchannel.rho_liquid[i,j]
        else:
            v1 = self.subchannel.vel_vapor[i,j]
            v2 = self.subchannel.vel_vapor[i,j+1]
            rho = self.subchannel.rho_vapor[i,j]
            
        # 计算湍流混合动量
        return beta * rho * abs(v1 - v2) * (v2 - v1)
        
    def _calculate_saturation_temperature(self, P):
        """根据压力计算饱和温度
        
        Args:
            P: 压力 (Pa)
            
        Returns:
            饱和温度 (K)
        """
        # 这里使用简化的饱和温度计算公式
        # 实际应用中应该使用更准确的热力学关系
        P_MPa = P / 1e6  # 转换为MPa
        return 273.15 + 100 * P_MPa / 0.101325  # 简化的线性关系
        
    def _calculate_saturation_enthalpy(self, P, phase):
        """根据压力计算饱和焓值
        
        Args:
            P: 压力 (Pa)
            phase: 相态 ('liquid' 或 'vapor')
            
        Returns:
            饱和焓值 (J/kg)
        """
        # 这使用简化的饱和焓计算公式
        # 实际应用中应该使用更准确的热力学关系
        T_sat = self._calculate_saturation_temperature(P)
        if phase == 'liquid':
            return 4200 * (T_sat - 273.15)  # 简化的液相焓值
        else:
            return 4200 * (T_sat - 273.15) + 2260e3  # 简化的气相焓值（加上汽化潜热）
        
    def _calculate_phase_change(self, k, i, j):
        """计算相变质量流量 (eq.10-11)
        
        Args:
            k: 相态 ('l' 或 'v')
            i, j: 网格位置
            
        Returns:
            相变质量流量 (kg/m³·s)
        """
        # 获取流体状态
        T_l = self.subchannel.temp_liquid[i,j]
        T_v = self.subchannel.temp_vapor[i,j]
        P = self.subchannel.pressure[i,j]
        
        # 计算饱和温度和焓值
        T_sat = self._calculate_saturation_temperature(P)
        h_l = self._calculate_liquid_enthalpy(T_l)
        h_v = self._calculate_vapor_enthalpy(P, T_v)
        h_l_sat = self._calculate_saturation_enthalpy(P, 'liquid')
        h_v_sat = self._calculate_saturation_enthalpy(P, 'vapor')
        
        # 计算各种状态下的相变率
        Gamma_l_sh = 0  # 液相过热
        Gamma_v_sh = 0  # 气相过热
        Gamma_l_sc = 0  # 液相过冷
        Gamma_v_sc = 0  # 气相过冷
        
        # 液相过热/过冷判断
        if T_l > T_sat:  # 液相过热
            C_l_sh = 0.1  # 液相过热传热系数
            Cp_l = self.subchannel.cp_liquid[i,j]
            Gamma_l_sh = C_l_sh * self.subchannel.rho_liquid[i,j] * Cp_l * (h_l - h_l_sat)
        elif T_l < T_sat:  # 液相过冷
            C_l_sc = 0.1  # 液相过冷传热系数
            Cp_l = self.subchannel.cp_liquid[i,j]
            Gamma_l_sc = C_l_sc * self.subchannel.rho_liquid[i,j] * Cp_l * (h_l_sat - h_l)
            
        # 气相过热/过冷判断
        if T_v > T_sat:  # 气相过热
            C_v_sh = 0.1  # 气相过热传热系数
            Cp_v = self.subchannel.cp_vapor[i,j]
            Gamma_v_sh = C_v_sh * self.subchannel.rho_vapor[i,j] * Cp_v * (h_v - h_v_sat)
        elif T_v < T_sat:  # 气相过冷
            C_v_sc = 0.1  # 气相过冷传热系数
            Cp_v = self.subchannel.cp_vapor[i,j]
            Gamma_v_sc = C_v_sc * self.subchannel.rho_vapor[i,j] * Cp_v * (h_v_sat - h_v)
            
        # 计算净相变率
        Gamma = Gamma_l_sh + Gamma_v_sh - Gamma_l_sc - Gamma_v_sc
        
        return Gamma if k == 'l' else -Gamma
        
    def _calculate_liquid_enthalpy(self, T):
        """计算液相焓值 (J/kg)
        
        Args:
            T: 温度 (K)
            
        Returns:
            液相焓值
        """
        # 以0℃液态水为基准
        Cp_l = 4200  # 假设恒定比热容
        return Cp_l * (T - 273.15)
        
    def _calculate_vapor_enthalpy(self, P, T):
        """计算气相焓值 (J/kg)
        
        Args:
            P: 压力 (Pa)
            T: 温度 (K)
            
        Returns:
            气相焓值
        """
        # 计算饱和温度
        T_sat = self._calculate_saturation_temperature(P)
        
        # 以0℃液态水为基准，加上汽化潜热和过热焓
        h_fg = self._calculate_latent_heat(T_sat)  # 汽化潜热
        h_f = self._calculate_liquid_enthalpy(T_sat)  # 饱和水焓值
        Cp_v = 2000  # 假设恒定比热容
        
        return h_f + h_fg + Cp_v * (T - T_sat)
        
    def _calculate_interfacial_heat_transfer(self, i, j):
        """计算界面传热
        
        Args:
            i, j: 网格位置
            
        Returns:
            界面传热量 (W)
        """
        # 获取温度
        T_l = self.subchannel.temp_liquid[i,j]
        T_v = self.subchannel.temp_vapor[i,j]
        T_sat = self._calculate_saturation_temperature(self.subchannel.pressure[i,j])
        
        # 计算界面面积
        alpha_l = 1 - self.subchannel.void_fraction[i,j]
        alpha_v = self.subchannel.void_fraction[i,j]
        d_h = config.GEOMETRY_PARAMS['hydraulic_diameter']
        A_i = 4 * alpha_l * alpha_v / d_h
        
        # 计算界面传热系数
        h_i = 1000  # W/m²K, 可以根据实际情况调整
        
        # 计算界面传热量
        q_i_l = h_i * A_i * (T_sat - T_l)  # 液相界面传热
        q_i_v = h_i * A_i * (T_sat - T_v)  # 气相界面传热
        
        return q_i_l + q_i_v
        
    def _calculate_heat_transfer_coefficient(self, k, i, j):
        """计算对流换热系数 (eq.6-9)
        
        Args:
            k: 态 ('l' 或 'v')
            i, j: 网格位置
            
        Returns:
            传热系数 (W/m²K)
        """
        # 获取温度
        T_w = self.subchannel.temp_wall[i,j]
        T_sat = self._calculate_saturation_temperature(self.subchannel.pressure[i,j])
        T_l = self.subchannel.temp_liquid[i,j]
        
        if k == 'l':  # 液相
            # 判断沸腾状态
            if T_w < T_sat:  # 单相液态
                # 使用Dittus-Boelter相关式
                Re = self.subchannel.rho_liquid[i,j] * abs(self.subchannel.vel_liquid[i,j]) * \
                     config.GEOMETRY_PARAMS['hydraulic_diameter'] / self.subchannel.mu_liquid[i,j]
                Pr = self.subchannel.cp_liquid[i,j] * self.subchannel.mu_liquid[i,j] / \
                     self.subchannel.k_liquid[i,j]
                     
                # 防止除零
                Re = max(Re, 1e-6)
                Pr = max(Pr, 1e-6)
                
                Nu = 0.023 * Re**0.8 * Pr**0.4
                h_l = Nu * self.subchannel.k_liquid[i,j] / config.GEOMETRY_PARAMS['hydraulic_diameter']
                
            elif T_l < T_sat:  # 过冷核态沸腾
                # 使用Thom相关式结合冷凝修正
                F_gam = 0.5  # 气泡冷凝系数
                P = self.subchannel.pressure[i,j]
                
                # 防止除零和数值不稳定
                dT = max(T_w - T_l, 1e-6)
                h_l = (1 - F_gam) * math.exp(2 * P / 1260 * (T_w - T_sat) / (2722 * dT))
                
            else:  # 饱和核态沸腾
                # 使用Thom相关式
                P = self.subchannel.pressure[i,j]
                dT = max(T_w - T_l, 1e-6)
                h_l = math.exp(2 * P / 1260 * (T_w - T_sat) / (2722 * dT))
                
            # 限制传热系数的范围
            h_l = max(min(h_l, 1e5), 1e2)  # 限制在100-100000 W/m²K之间
            return h_l
            
        else:  # 气相
            if self.subchannel.void_fraction[i,j] < 0.999:  # 非单相气态
                return 0  # 忽略壁面与气泡之间的传热
            else:  # 单相气态
                # 使用Dittus-Boelter相关式
                Re = self.subchannel.rho_vapor[i,j] * abs(self.subchannel.vel_vapor[i,j]) * \
                     config.GEOMETRY_PARAMS['hydraulic_diameter'] / self.subchannel.mu_vapor[i,j]
                Pr = self.subchannel.cp_vapor[i,j] * self.subchannel.mu_vapor[i,j] / \
                     self.subchannel.k_vapor[i,j]
                     
                # 防止除零
                Re = max(Re, 1e-6)
                Pr = max(Pr, 1e-6)
                
                Nu = 0.023 * Re**0.8 * Pr**0.4
                h_v = Nu * self.subchannel.k_vapor[i,j] / config.GEOMETRY_PARAMS['hydraulic_diameter']
                
                # 限制传热系数的范围
                h_v = max(min(h_v, 1e4), 10)  # 限制在10-10000 W/m²K之间
                return h_v